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Abstract 

In this paper we describe a novel approach for the solution of inviscid flow problems for subsonic 
compressible flows. The approach is based on canonical forms of the equations, in which subsystems 
governed by hyperbolic operators are separated from those governed by elliptic ones. The discretizations 
used as well as the iterative techniques for the different subsystems, are inherently different. Hyperbolic 
parts, which describe, in general, propagation phenomena, are discretized using upwind schemes and 
are solved by marching techniques. Elliptic parts, which are directionally unbiased, are discretized 
using h-elliptic central discretizations, and are solved by pointwise relaxations together with coarse grid 
acceleration. The resulting discretization schemes introduce artificial viscosity only for the hyperbolic 
parts of the system; thus a smaller total artificial viscosity is used, while the multigrid solvers used are 
much more efficient. Solutions of the subsonic compressible Euler equations are achieved at the same 
efficiency as the full potential equation. 
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1 Introduction 


In the past decade a substantial effort has been invested in solving the Euler equations, with multigrid 
methods playing a central role. The two major directions of research in multigrid solution of the Euler 
equations are the use of coarse grids to accelerate the convergence of the fine grid relaxations, and the use 
of defect correction as an outer iteration, coupled with use of multigrid to solve for the low order operator 
involved [3]. Extensive research has been conducted in both directions. Unfortunately, both approaches 
can be shown to have limited potential. Methods based on defect correction have /i-dependent convergence 
rates for hyperbolic equations, where h is the mesh spacing; other multigrid methods have p-dependent 
convergence rates, where p is the order of the scheme involved. This unacceptable situation motivated the 
research outlined here. 

The poor behavior of coarse grid acceleration for hyperbolic equations, which becomes even worse with 
high order discretizations, leads us to conclude that coarse grids should not be used to accelerate convergence 
for hyperbolic problems. Rather, the relaxation process alone should converge all components of such 
problems. This is possible, since hyperbolic problems describe propagation phenomena for which marching 
in the appropriate direction is a highly effective solver. For elliptic problems, on the other hand, local 
relaxation with good smoothing properties can be greatly accelerated by coarse grid correction. Moreover, 
elliptic problems cannot be solved efficiently by local relaxation alone, so that coarse grid acceleration is 
essential. In short, hyperbolic equations do not need coarse grid acceleration, while elliptic equations require 
such acceleration. 

These observations have motivated a study concerning the separation of the hyperbolic and elliptic parts 
in steady state inviscid flow calculations. The result is a canonical form for the inviscid equations, where the 
hyperbolic and elliptic parts reside in different blocks of an upper triangular matrix for the system [5]. These 
canonical forms are the analog of the decomposition of the time dependent one-dimensional Euler equations 
into characteristic directions and Riemann invariants. The insight gained by the use of the canonical variables 
enables one to construct genuinely multidimensional schemes for these equations. It unifies the treatment of 
the compressible subsonic case with the incompressible case, although these two cases have been studied by 
different methods until now. Canonical boundary conditions [5] which enable the proper numerical treatment 
of general boundary conditions are also obtained. Moreover, these canonical forms enable the construction 
of new solvers having much better efficiency than existing solvers. 

The new discretization schemes, which are based on these canonical forms, use upwind discretization 
only for the hyperbolic, variables, and use a central h-elhpttc discretization for the elliptic variables. This 
gives a better representation for the physical phenomena, since elliptic problems do not have a bias in 
any spatial direction, a property that should hold true for the discretization as well, if possible. Upwind 
discretization for hyperbolic problems, on the other hand, is compatible with the bias of information flow in 
the physical problem. The canonical form, therefore, allows for a better treatment of the inviscid equations. 
The resulting schemes are also compatible with the uniqueness properties of the inviscid equations under 
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different geometries and boundary conditions. In particular, the nonuniqueness of solutions for exterior flows 
around smooth bodies is clearly evident with these schemes; only the addition of a global condition, such 
as circulation at infinity, insures uniqueness. In existing schemes this issue is much more obscure and there 
seems to be no direct analog of the physical behavior. 

We discretize the elliptic part of the system with staggered grid schemes, which are h- elliptic. The 
importance of h-elliptic discretizations is that they yield solutions free of spurious or weakly spurious modes, 
in the limit as the Mach number goes to zero. The construction of fast solvers for such discretizations is 
well understood, even in the case of systems, and is quite straightforward. By contrast, the non-staggered 
schemes commonly used would require special techniques for the construction of optimally efficient multigrid 
solvers (having efficiency comparable to that of the full potential equation). Such solvers require multiple 
coarse grids on each level, even for regular grids having cell aspect ratios near one [2]. 

The elliptic and hyperbolic parts of the equations are treated differently by the solver. Unlike existing 
solvers, which use coarse grids to accelerate the hyperbolic part of the system as well, the new method 
computes the hyperbolic part via relaxations based on the canonical forms. This involves marching in the 
streamwise direction for the hyperbolic quantities, which are the entropy s and total enthalpy H . The rest 
of the unknowns, e.g., the velocity components, are relaxed by a Kaczmarz relaxation using preconditioned 
residuals, yielding a smoothing rate identical to that for the full potential equation. 

Numerical results are given for a two dimensional flow around an ellipse and flow in nozzle. These 
problems already include the major difficulties in real problems and serve as a good test for the method 
proposed. Second order schemes are used for both cases and the solutions are obtained with convergence 
rates similar to that of the full potential equation (although the work involved here is larger accounting for 
the multiple equations). 

2 Canonical Forms and Discretization Rules 

The discretization and efficient solution of elliptic systems of partial differential equations is quite well 
understood. One of the important concepts here is h-ellipticity [1]. It guarantees that the stability of high 
frequencies for the discrete problem is in correspondence to that of the differential system. For the latter, 
ellipticity is defined in terms of the symbol P(u>) as 

detP(w)>CV| 2m , (2.1) 

while h-ellipticity is defined as 

det P h {9) > C\0\ 2m , \0\ < 7T (2.2) 

Discretizations which are h-elliptic admit local relaxation methods with good smoothing properties. This, 
together with efficient coarse grid acceleration for smooth components, makes standard multigrid methods 
very efficient for such discretizations. Although other types of discretization also admit fast multigrid solvers, 
we restrict our focus to h-elliptic discretization for elliptic problems. 
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The discretization of hyperbolic equations, which usually describe propagation phenomena, can be done 
naturally using upwind biased schemes. However the application of the above ideas to the steady state 
inviscid incompressible and subsonic compressible equations is not straightforward, since these equations are 
neither elliptic nor hyperbolic, but rather mixed hyperbolic-elliptic. The optimal treatment of the problem 
should therefore include an identification of these two parts, which have inherently different behavior and 
call for different numerical treatment both on the level of the discretization and the solver. The device for 
this is a canonical form of the equations, described in detail in [5]. 

For a two dimensional flow, the canonical form of the compressible Euler equations is 

/ D\ Do 0 0 

qD y —qD x — y( y _ i ) Do 
“1 6 -TpQ 0 

V 0 0 0 pQ 

where: 

D\ = p/c 2 ((c 2 — u 2 )D x — uvD y ) 

D 2 = p/c 2 {(c 2 - v 2 )D y - uvD x ) (2.4) 

Do — vD x xiD y 

In view of these canonical forms we can use the following discretization rule for the invisicd equations: 

1. Use central (unbiased) h-elliptic. discretizations for elliptic subsystems 

2. Use upwind biased schemes for hyperbolic subsystems 

3 Discretization 

Let a domain ft £ JR 2 be divided into arbitrary cells. Let the vertices, edges, and cells be V,E, and C 
respectively. The well known Euler formula 

#V + #C + #Ao/e* = P+1 (3.1) 

suggests several possibilities for discretization of different systems on structured and unstructured meshes. 
Rewriting the Euler formula as 

#V r + #V + #c + #C + holes = + 1 (3.2) 

one obtains the following choice of discretization. Let H be associated with the cell centers, the normal 
velocity components be at the edges, and the entropy be at the vertices. Quantities other than these 
are calculated by well known algebraic relations for the thermodynamic quantities and by averaging for 
the tangential velocity components. With each cell, we associate one continuity equation and one energy 
equation, while the two momentum equations are associated with each vertex. The following diagram is then 
obtained: 
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Figure 1: 
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(3-3) 


The control volumes for an unstructured mesh, for the continuity, energy, and momentum equations are 
shown in Figure 2. 

The canonical form for the compressible equations suggests that only the entropy and the total enthalpy 
will be discretized using upwind biased schemes, and only in the appropriate terms; that is, only those terms in 
which a derivative in the streamwise direction occurs will be discretized with upwind bias. Other derivatives 


of these quantities will be discretized using central differencing. Let ei = (u/q y v/q) ) and ei = (— v/q.ujq). 
Decomposing the unit normal vectors to the edges as 


n = (n ei)ei -f (n • e 2 )e 2 (3.4) 

we obtain the following discretization of the pressure terms 


d P = X^ P ( n < ' e i ) e i + Pi( n > ‘ e 2) e 2l ds i (3.5) 

(3.6) 
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where 



P up = p(S up , H c , q 2 ), 

(3.7) 


p c =p(S e ,H e ,q 2 ). 

(3.8) 

and superscript 

refers to an upwind biased approximations, while superscript ” c” 

refers to a central 

approximation. 

The full discretization is then 



^ ^ Pj ^ j ‘ n j j = 0 > 
j €7 

(3.9) 


•n,)V l rf5, + dp=0, 

/€ 7 

(3.10) 


=0, 

(3.11) 


where p is computed using a symmetric formula at subsonic points, and using an upwind biased one for 
supersonic points. All other quantities involved use central approximations. Control volumes for the different 
equations are denoted by 7 and 7 . Thus the only quantities that have been upwinded are the entropy and 
total enthalpy. This results in a reduced artificial viscosity. 

4 Multigrid Algorithm 

The multigrid solver, like the discretization, is based on the canonical forms mentioned in section 2. The 
main ingredient differing from other multigrid methods is the relaxation scheme. Other elements of the 
multigrid method are standard and will be mentioned only briefly. 

As can be seen from the canonical form, the hyperbolic and elliptic parts for subsonic flows are separated. 
Since these subsystems are of very different nature, it is unlikely that the same numerical process will be 
optimal for both. Indeed, it can be shown that coarse grids are inefficient in accelerating certain smooth 
components for hyperbolic problems. For these components convergence rate is roughly (2 P — 1)/2 P for a 
p — order method. The better the scheme, the less coarse grids help! 

This behavior suggests that coarse grids are not appropriate for accelerating convergence for hyperbolic 
problems. The relaxation should therefore converge all hyperbolic quantities in the problem. While for 
elliptic problems relaxation cannot be efficient for smooth components, for hyperbolic problems the situation 
is different. Marching in the direction of the physical flow is very efficient in converging all hyperbolic 
components eliminating the need for coarse grid acceleration. For elliptic problems, on the other hand, 
relaxation techniques with good smoothing properties, combined with coarse grid acceleration, yield optimal 
solvers. The separation of the different subsystems presented by the canonical form allows one to construct 
an optimal solver for the full system. Marching techniques will be used for the hyperbolic quantities, while 
local relaxation with good smoothing will be used for the elliptic parts. 
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4.1 Relaxation 


Let the residual of the compressible Euler equations be denoted by (R p , R pu , R pv , R H ) and the ones for the 
canonical form by (r* , r 2 f r^ y r*). Then the following relation prevails 
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The relaxation for the Euler equations is employed as follows. The total enthalpy is relaxed first, using 
the preconditioned residual r* 


H y +—H y + (rt) y ./c“ (4.2) 

Next the entropy is relaxed using 

s y * Sy + {r^)yfcy } (4 3 ) 

where cif and Cy are the diagonal coefficient in the discrete version of uD x -f vD y and —pT(uD x + t>D y ), 
respectively. This is followed by relaxing the continuity and vorticity equations. The relaxation for the 
continuity equation is done by keeping the values for the density frozen, i.e., 

v> • n, ♦ — V, • n, + {r\) y pjdsj/ ^(ptds t ) 2 , j € 7, (4.4) 

and the vorticity equation is relaxed as: 

V|-t|«— V,-t t + tf/iqp^dst/Ylidsk) 2 lei (4.5) 

key 

Note that the preconditioning of the discrete system does not result in an exact upper triangular form. 
Lower diagonal terms exist, but these are of order 0(h 2 ) and do not affect the design of the relaxation and 
other numerical processes. 

The coarsening part of the multigrid method used is standard, so its details are omitted. Coarse grids are 
created by combining neighboring fine grid cells into a coarse grid cell. Linear interpolation of corrections 
and full weighting of residual and functions are used in an FMG-FAS formulation [1], 

5 Numerical Results 

We present here numerical results for subsonic, cases of flow in a nozzle and around an ellipse. These cases 
already present most of the difficulties encountered in subsonic compressible flows. For both cases body 

fitted grid were used. In the case of the ellipse the grid extended to to a distance of 7 chords. 

Figure 3 and 4 show results for flow in a nozzle. Mach lines as well as relative errors in entropy are shown. 
Cases of inflow Mach number of .2 and .02 are shown. Observe the small size of the errors in entropy, as 
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well as the symmetry of the solution. The convergence rate of the multigrid method used is shown in Figure 
5, in the top two pictures. A reduction of eight orders of magnitude in 10 cycles is obtained. Note that the 
rate of convergence is independent of the Mach number. 

Figure 5 shows results for a flow around an ellipse at a 10 degrees angle of attack, with free stream 
Mach number of 0.1, with zero circulation at infinity. This is a case of special difficulty. The exact values 
for the lift and drag coefficient are zero. The calculated values are Cl = 2.3 x 10“ 3 ,C/> = —7.7 x 10 3 . 
A grid refinement study shows that the lift and drag coefficient approach zero with mesh refinement. The 
convergence history for this problem is shown in the bottom picture of figure 6. The full FMG history 
(including coarse levels) is shown. The convergence rate in this case is slower than for the nozzle flow cases 
but still very fast compared with existing methods. 

6 Conclusion 

A new approach for the discretization and the solution of the Euler equations have been presented. It is 
based on a canonical form of the equations which allows for discretization which involve upwind biased 
discretization only for the physically biased quantities, that is, entropy and total enthalpy. The multigrid 
method used for that disctetization shows essentially optimal convergence rates for the second order schemes 
used. Moreover, nonlifting solutions around smooth bodies at angle of attack are obtained. Unlike most 
other methods, the performance of the method does not degrade as the Mach number approaches zero. 
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Figure 5: Residulas History 
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